# Functions to run caliper test

caliper_test <- function(tstats, caliper, threshold = qnorm(.975)) {
  n_window <- sum(abs(tstats) >= threshold * (1 - caliper) &
                    abs(tstats) <= threshold * (1 + caliper))
  n_above <- sum(abs(tstats) >= threshold &
                   abs(tstats) <= threshold * (1 + caliper))
  data.frame(
    caliper = caliper, below = n_window - n_above, above = n_above,
    pvalue = pbinom(n_above - 1L, n_window, 0.5, lower.tail = FALSE)
  )
}

run_all_caliper_tests <- function(tstats) {
  do.call(rbind, lapply(c(0.1, 0.15, 0.2), caliper_test, tstats = tstats))
}

format_caliper_test_table <- function(caliper_results) {
  caliper_results$caliper <- paste0(100 * caliper_results$caliper, "\\%")
  caliper_results$pvalue <- format(caliper_results$pvalue, digits = 2)
  outtab <- as.matrix(caliper_results)
  outtab <- format(rbind(c("Caliper", "Below", "Above", "$p$-value"), outtab))
  outtab <- apply(outtab, 1, paste0, collapse = " & ")
  outtab <- paste0(outtab, " \\\\")
  outtab <- c(
    "\\toprule",
    outtab[1],
    "\\midrule",
    outtab[2:4],
    "\\bottomrule"
  )
  paste0(outtab, collapse = "\n")
}

make_caliper_test_table <- function(tstats) {
  format_caliper_test_table(run_all_caliper_tests(tstats))
}

# Table 1: Caliper test with all main t-stats

rd_papers <- read.csv("data/cl-rd-papers-meta.csv")
rd_papers <- subset(rd_papers, main_estimate)

caliper_test_table <- make_caliper_test_table(rd_papers$tstat)
cat(caliper_test_table, file = "output/tables/table1-main-caliper-test.tex")


# Table S4: Caliper tests by BW selection method

rd_papers_autoBW <- subset(rd_papers, bw_select_strat != "arbitrary")
rd_papers_not_autoBW <- subset(rd_papers, bw_select_strat == "arbitrary")
stopifnot(nrow(rd_papers_autoBW) + nrow(rd_papers_not_autoBW) == nrow(rd_papers))

caliper_test_table <- make_caliper_test_table(rd_papers_autoBW$tstat)
cat(caliper_test_table, file = "output/tables/tableS4a-auto-caliper-test.tex")

caliper_test_table <- make_caliper_test_table(rd_papers_not_autoBW$tstat)
cat(caliper_test_table, file = "output/tables/tableS4b-noauto-caliper-test.tex")


# Table S5: Caliper tests for re-analysis t-stats

# Load reanalysis data
rd_papers <- read.csv("data/cl-rd-papers-meta.csv")
reanalysis_res <- read.csv("data/reanalysis-results.csv")

col_to_merge <- c("estimate_code", "bw_select_strat", "main_estimate")
reanalysis_res <- merge(
  reanalysis_res,
  subset(rd_papers, data_avail, col_to_merge),
  by = "estimate_code"
)
rm(rd_papers, col_to_merge)

# Subset to main estimates and those not using CCT
reanalysis_res <- subset(reanalysis_res, main_estimate & bw_select_strat != "cct mse optimal")

caliper_test_table <- make_caliper_test_table(reanalysis_res$cct_conven_tstat)
cat(caliper_test_table, file = "output/tables/tableS5a-cct-conven-caliper-test.tex")

caliper_test_table <- make_caliper_test_table(reanalysis_res$cct_robust_tstat)
cat(caliper_test_table, file = "output/tables/tableS5b-cct-robust-caliper-test.tex")
